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We study the ground state and the first three radially excited states of a self-gravitating Bose-Einstein- 
Condensate with respect to the influence of two external parameters, the total mass and the strength of in¬ 
teractions between particles. For this we use the so-called Gross-Pitaevskii-Newton system. In this context we 
especially determine the case of very high total masses where the ground state solutions of the Gross-Pitaevskii- 
Newton system can be approximated with the Thomas-Fermi limit. Furthermore, stability properties of the 
computed radially excited states are examined by applying arguments of the catastrophe theory. 


I. INTRODUCTION 


The first Bose-Einstein-Condensate (BEC) was produced in 
1995 by the group of E. Cornell and C. Wiemann m , 70 years 
after its prediction by S. Bose and A. Einstein 1925 El- Be¬ 
ing a realization of a macroscopic quantum object, BECs have 
several interesting properties and therefore have been exten¬ 
sively studied experimentally as well as theoretically. 

BECs are a dilute quantum gas with short range dipole in¬ 
teractions between the atoms. Thus, for ultracold temper¬ 
atures they are described by means of the Gross-Pitaevskii 
equation (GP equation) Em 

o 12 

( 1 ) 

coupled to an external potential V describing in particular the 
traps (e.g. a harmonic potential) El- Here ip represents the 
wave functions of the condensate. The parameter g describes 
the self-interaction. Despite their diluteness, it is interesting to 
discuss self-gravitating BECs from a conceptual point of view 
and also from an experimental and astrophysical perspective. 

Self-gravitating quantum systems have been proposed by 
R. Penrose in his discussion of quantum state reduction by 
gravity ||6l. He considered a self-gravitating Schrodinger field 
described by the Schrodinger-Newton (SN) equations 

d fi^ 

ih—Ip = - — Alp + V'ip + m^'ip 

ot 2m (2) 

A$ = 47rG|V'P, 

where G is the Newton constant. This setting would dramat¬ 
ically change the concept of quantum mechanics, where one 
only assumes interactions with other external helds. 

In an astrophysical context, R. Ruffini and S. Bonazzola 
were the first to discuss self-gravitating bosons which are ex¬ 
clusively trapped in their own gravitational potentials, as a 
concept for boson stars 0- The SN equations have been stud¬ 
ied further extensively by R. Harrison and I.M. Moroz et al. 
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Here we are going to discuss self-gravitating BECs given 
by the Gross-Pitaevskii-Newton (GPN) system 


d h? 

=- 7 ;—^'<P+ 'rn^'ip+ Vtp + g\'ip\^'ip 
at 2m (3) 

A$ = 47rG|V'P- 


Contrary to the SN equations, a self-gravitating BEC will 
still have a meaning within the standard quantum mechani¬ 
cal framework. This is as well the case for the limit p —^ 0, 
which recovers the SN equations in 0. In this case however 
^ in 0 has to be interpreted as wave functions of the con¬ 
densate. On the experimental side, one may think of creating 
a self-gravitating BEC through roque waves which can be ex¬ 
cited within BEC due to the nonlinearity of the GP equation 
inoiini. The GPN system also may serve as model for nonrel- 
ativistic boson stars. The SN equations are the non-relativistic 
limit of the Einstein-Klein-Gordon (EKG) system ii. Anal¬ 
ogously, the GPN system can be obtained as nonrelativistic 
limit of a generalized EKG system 


R 


Dtp + U{tp) = 0 


(4) 


with an extra potential U for the Boson field describing self¬ 
interactions due to local interactions between the Bosons. 

is the energy momentum tensor of the Klein-Gordon 
field. Usually U{tp) is given by some polynomial. Eor the 
usual Klein-Gordon equation without local self-interactions 
we have U(tp) = mtp. Such EKG systems have been exten¬ 
sively studied as model for relativistic boson stars, see oini 
for reviews. 

Giant self-gravitating BECs have been suggested as candi¬ 
dates for dark matter (DM) halos (e.g. niH^). While C. G. 
Bohmer and T. Harko discussed the Thomas-Eermi limit (TE 
limit), P.-H. Chavanis discussed the full GPN system numeri¬ 
cally ED and analytically EOll for ground state solutions. 

The nature of DM is one of the major quests in cosmol¬ 
ogy and theoretical physics. Among others the flatness of ob¬ 
served galaxy rotation curves can not be explained by New¬ 
tonian gravity or standard general relativity, if only the visi¬ 
ble matter is considered. As a result one assumes the exis¬ 
tence of invisible dark matter which forms a spherical halo 
around galaxies. There are a number of suggested dark mat¬ 
ter models, most popular the A cold dark matter (ACDM) 
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model which comprises weakly interacting massive particles 
(WlMPs). But there arise difficulties to explain the observed 
distribution of the invisible matter at galactic centers (in the 
scales of the order of 1 kps and smaller) with the CDM model. 
It leads to cuspy density profiles im instead. In the context 
of the CDM model self-gravitating BECs are discussed as a 
dark matter candidate which solves the occurring cusp prob¬ 
lem. The quantum properties of self-gravitating BECs lead 
to a repulsive force due to the Heisenberg uncertainty princi¬ 
ple which prevents the forming of cusps in the density profile 
II 22 I . Rotation curves induced by self-gravitating BECs have 
been compared to observed galaxy rotation curves using dif¬ 
ferent relativistic ini and non-relativistic models (e.g. using 
the TE limit ifTSl . neglecting particle interactions ifThl or dis¬ 
cussing rotating self-gravitating BECs ll23l ) and seem to ap¬ 
proach the observed curves. 

Now the question arises, whether the obtained solutions 
for a self-gravitating BEC are actually stable. The exis¬ 
tence of a sort of Jeans instability for an infinite spatially 
homogeneous distribution of self-gravitating bosons could be 
shown in Esiia and in (taking into account short- 
range interactions). Discussing self-gravitating BECs, it was 
shown in case of a relativistic treatment that for a poten¬ 
tial [/('!/') of the form Ui^tp) = X — aip'^ + or 
U2i'4’) = A'('0'* — a''0^ + for A, A', a, a', 6,6' > 0, 
stable and unstable ground state solutions exist 11611271. If 
self-interactions between particles are neglected the found 
ground state solutions are stable |28l. Radially excited states 
on the other hand seem to be unstable either way, if the po¬ 
tential U 2 or no self-interaction is considered ll28l |2^ . In the 
non-relativistic limit, the SN equations and the GPN system 
are discussed. Ground state solutions of the SN equations are 
shown to be stable, while radially exited states of the system 
were found to be unstable ii. In case of the GPN system, sta¬ 
ble and unstable ground state solutions exist EOl [3011311 . but 
the radially excited states seem to be unstable ll32l . However 
the investigation of excited states is still of interest. In case 
of the SN equations S.-J. Sin had to use excited states with 
more than four nodes to get a satisfying approach of the ob¬ 
served rotation curves. Eurthermore, L.A. Urena-Lopez and 
A. Bernal found that a superposition of the ground state and 
excited states can be stable 13^ . 


In this paper we continue the work of P.-H. Chavanis 
by calculating and analyzing the first three radially excited 
states (wave function solutions with a characteristic number 
of nodes) of the GPN system for influence of the two external 
parameters: total mass and strength of the particle interaction. 
We restrict to non-rotating systems. 


This paper is organized as follows. In section [I^ the GPN 
system is introduced as well as two of its limits: the TP limit 
and the non-interacting limit. In section III the used numer¬ 
ical procedure is explained and confirmed to work properly 
by comparing computed solutions with solutions published by 
other authors. In section IV the computed results of the first 
three radially excited states are compared to the numerical re¬ 
sults for the ground state. Ground state solutions were already 
discussed by P.-H. Chavanis ||T9[ . but are nevertheless com¬ 
puted in this paper as well for reasons of completeness. In 


section |V] the stability properties of the found solutions are 
discussed. As already mentioned above, the radially excited 
states of the GPN system are found to be intrinsically unsta¬ 
ble by P. S. Guzman and L. A. Urena-Lopez by studying 
the time evolution of equilibrium solutions, while allowing a 
flow of particles out of the numerical domain. However the 
stability of the GPN system was discussed only for a few dif¬ 
ferent values of the parameter g, which indicates the strength 
of particle interactions. In this paper we use a different way 
to approach the stability question. By applying arguments of 
the catastrophe theory we are able to give statements about a 
whole branch of possible configurations for a self-gravitating 
BEC. Pinally, in section |VI| we draw some comments and 
conclusions. 


II. THE GROSS-PITAEVSKII-NEWTON SYSTEM 


A BEC in its own gravitational held is described by the 
GPN system lHHIlOl 

ot 2m (5) 

-f 5 |T'(f,f)|^T'(r,f) 

A<l>(r,f)= AttCto f)|^ . (6) 

With the self-interaction factor g — dTrS^a/m short range 
interactions are taken into account. Here a is the scattering 
length of the particles forming the BEC which can be chosen 
either positive (repulsive) or negative (attractive), and m is the 
mass of each single particle which is part of the BEC. 

The wave function Ar (r, t) satisfies the normalization con¬ 
dition 


N = J d^r \^{r,t)\\ (7) 

where N is the total number of particles. |T'(r, f)|^ corre¬ 
sponds to the particle density and p{f) = m\'i’{r,t)\ gives 
the mass density. Thus, the total mass M of the BEC can be 
computed with 


M = mN = m 


d^f |4'(r, t)|^ 


( 8 ) 


We obtain the time-independent GPN system for stationary 
solutions of the form T'(r, t) = ip{f) ex-p{—iEt/h), 


Eij;(r) = - V'^'ipir) + m^(r)tfj{f) + g |0(r)|^ iplE) 

2m 

(9) 

A$(f^ = dTrGm |'0(f*)l^ ) (10) 

with E being the eigenenergy of the GPN system. 
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A. The energy functional 

Now we introduce the funcional i?tot[V'] which is associated 
with the total energy of the GPN system; 

" ( 11 ) 

+ 1 

Here <i) is a function according to An extremum 

of the total energy at fixed total mass M is given by the varia¬ 
tional principle 



— neglected. This limit was studied 

in iflSl . Due to the dominant resulting repulsive inter¬ 
actions (solutions in the TF-limit can only be found for 
p > 0) the self-gravitating EEC is prevented from grav¬ 
itational collapsing. In the TF limit we obtain an ap¬ 
proximate solution for the ground state of the GPN sys¬ 
tem. This is due to the fact that in the TF limit solutions 
can only be found if there exists a repulsive force F'si(fO 
for all r. This force is caused by the self-interaction 
term p |'0(r) is given by 

= -Vp{r). (16) 

m 


SEtot — aS J cPf |T'(r,f)|^ — Mj = SI = 0 , (12) 

where equation ([^ is used, a is a Fagrange multiplier which 
takes into account the mass constraint and which can be inter¬ 
preted as a chemical potential. The functional 


A repulsive force Fsi{r) < OVr can be found only for 
the ground state for which the condition V p{r) < 0 Vr 
is satisfied. Excited states do not fulhll this condition. 

The density prohle p{r) of a spherical symmetric GPN 
system in the TF limit is given by lIMIl 


I = J IV'(^r 




(13) 


is called the energy functional and its variation 51 is given by 


p(r) = sin f , for r < Rtp , (17) 

7rr \RtfJ 

where po is the central density at r = 0, and Rjp = 
TT{ah‘^ is the radius of the EEC where p{r) 

becomes zero. The eigenenergy E is determined as 


51 = J — (A'0i5'0* -f At/)*(5-0) 

+ (m<i) -I- g — a) + 'tp*5'ip) . 


(14) 


Since $ depends on ip, one has to variate $ in ( [T3]) as well. 
This results in the loss of a 1/2 factor from (0 to (|l41i at the 
<i)-term. Ey identifying a as eigenenergy E, solving 51/Sip = 
0 or 51 /Sip* = 0 results in the time-independent GPN system 
s and ( fT0| ). Wave functions, which correspond to extrema 
of the energy functional are therefore solutions to the GPN 
system. 


GMm 

Rtp 


(18) 


III. SOLVING THE GPN SYSTEM 


As we are interested in radially exited states only, we 
choose a spherically symmetric ansatz for the wave function, 
ip{r) = f{r). In the following we use E = huj. 


A. Rescaling 


B. Limiting cases of the GPN system 

We now describe two threshold regions of the solution man¬ 
ifold of the GPN system which we will need later. 

(i) In the non-interacting limit the self-interaction term in 
can be neglected 


For solving the GPN system it is convenient to rescale the 
variables and make them dimensionless. The introduction of 
natural length and energy units f = r and w = ^^"^5 w 
results in the rescaling of the other variables. With a further 
dimensionless scaling factor A with A G {K | A > 0} we ob¬ 
tain 


5|0(f)|"0(f) ^0. (15) 

Then the GPN system reduces to the SN system. The 
non-interacting limit is obtained for p = 0 or very small 
total masses M (i.e. |'0(r)|^ —)■ 0). 

(ii) The Thomas-Fermi limit (TF limit) is characterized by a 
large number of particles so that the kinetic term in (|^, 


9 = 

’"2m4GA2^’ 

f{r) 

UJ = 


$(r) 


fF ^ 

M 


2rrFGX'' 


1 


\ IF J 

A2$(r), 


2G^rrF 


fF 


( 19 ) 
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This rescaling results in the following dimensionless time- 
independent GPN system 


^ = kr)hr) - w/(f) 

dr^ 

<P^{f) ^ - 2 _ 2 #(f) 

df"^ f df 


2dm 

f df 



The scaling factor A does not appear in ( |20| ) showing that A 
leaves ( |20| ) form invariant. This characteristic feature can be 
used to obtain solutions for different values of g and M di¬ 
rectly from already computed solutions, by changing the value 
of A. 

In the following, the dimensionless GPN system ( |20l i will 
be solved numerically by using the FORTRAN program COL- 
SYS. 


B. The numerical procedure 


COLSYS can be used for numerically solving mixed-order 
problems of systems of ordinary differential equation (ODE) 
with given boundary conditions. (More information can be 
found in ll^ .) In order to enable the computation of the 
eigenenergy oj of the dimensionless GPN system ( |20l i an ad¬ 
ditional differential equation has to be added 



( 21 ) 


This equation completes the set of ODEs. 

The boundary conditions for the system of ODEs are 
chosen as follows: 


a) /'(0)=0, d) /(i?)=0, 

b) ^'(0)=0, e) = (22) 

M 


c) 6j\R)= 0, f) /(O) = /, or </.'(!?) = 


i?2 ’ 


where the prime denotes d/df. R represents a fixed and very 
large value of f which will be determined later. It is used to 
fix the behavior of the solution far away from the origin. The 
boundary condition f) can be used to preset either a value for 
/ (0) or the total mass M of the EEC. 


We start our computations with solutions (ground state and 
the first three radial excitations) characterized by p = 0 and 
M = 1. The obtained solutions are compared to those com¬ 
puted in m (in fact, by our method we were able to repro¬ 
duce the solutions). Branches of ground state solutions and 
branches of solutions of the first three radially excited states, 
respectively, can be computed by slowly changing the value 
of g or of M. A full analysis of the influence of the param¬ 
eters g and M on the GPN system needs to either set a fixed 
value of M and vary g or to vary M for a fixed g > 0 and a 
fixed g < 0. Solutions for other combinations of {g, M) are 
then obtained from these computed solutions by the variation 
of the value of the scaling factor A. 

Eor a comparison of the obtained solutions we match the 
values of the computed eigenenergy ui = E/H with the appro¬ 
priate values ujh calculated in 18]. Since in |8l N, H, G, m, A 
are set equal to 1, ujh has to coincide with 2uj (see scaling 
condition for w, equation ([Tg])). Table [^summarizes the com¬ 
parison of our results with those in IS]. Deviations are at most 
of order of 10“®. Thus, our ansatz and software seems to work 
correctly. Eurthermore the shapes of the computed ground 
state wave function and the corresponding gravitational po¬ 
tential for 5 = 0 and M = 1 are compared with the results 
obtained by R. Ruffini et al. in 121. The calculated functions 
perfectly coincide what further supports our conclusion. 


IV. SOLUTIONS 


Typical profiles of the obtained wave functions and their 
corresponding gravitational potentials of the GPN system are 
shown exemplarily in Pig. The number of nodes of a wave 
function is denoted k, the wave function itself is denoted fk 
and the corresponding gravitational potential is denoted (j)k. 
Eor the ground state we have k = 0 and the first three radially 
excited states are indicated hy k = 1,2,3. The solutions in 
Pig. [^correspond to g = M = 1. We observe that for larger k 
the maximum of fk{f) at the origin decreases, and that fk{f) 
more slowly approaches the abscissa. Thus, radially excited 
solutions of self-gravitating BECs have lower density in the 
center and a larger extension. Nodes of the radially excited 
states create “bulges” in their associated gravitational poten¬ 
tials. 


State 

U) 

2 ■ u) 

UJH in (U 


ground state 

-0.081384604 

-0.162769208 

-0.162769291 

5 • 10"’’ 

1. ex. state 

-0.015398269 

-0.030796537 

-0.030796561 

o 

1 

2. ex. state 

-0.006263051 

-0.012526101 

-0.012526108 

6 • 10"’’ 

3. ex. state 

-0.003373660 

-0.006747320 

-0.006747330 

1•10“® 


TABLE I. Computed values of the eigenenergy 6j of the ground state 
(fc = 0) and the first three radially excited states (k — 1, 2, 3) for 
g = 0, M — 1 compared with the values of the eigenenergy ujh 
computed in (8). Deviations A 2 a,uijj of 2a) from ujh do not exceed 
10 "®. 


A. The total mass - eigenenergy relation 

In Pig. [^ the eigenenergy Cjk of a wavefunction solution 
with k nodes is plotted as a function of the total mass M (a) 
for a fixed value g > 0 and (b) for a fixed value g < 0, re¬ 
spectively. These two plots are combined in a log-log plot in 
Pig. Eor g > 0 and g < 0 the behavior of ujk{M) shows 
the same qualitative shape for the ground state (fc = 0) and 
the first three radially excited states (fc = 1,2,3). Eor small 
masses M and g ^ 0 the function d)fc(M) looks like uj{M) in 
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FIG. 1. (a) Computed solutions for the wave function and (b) com¬ 
puted solutions for the gravitational potential of the ground state 
(fc = 0) and the first three radially excited states (k = 1, 2, 3) for 
g = M = l. 


the non-interacting limit g = 0. This behavior is described by 

(23) 


^k,c 

UJk = ^r-M 

M2 


6jk,c is the computed eigenenergy of a wave function solution 
with n nodes for M = M^- M^ can be chosen arbitralily. The 
relatio n ([2^ can be deduced from the rescaling relation of M 
and g (|19|l. 

For a fixed value g > 0, uJo{M) of the set of ground state 
solutions closely approaches the course of Cj{M) in the TF 
limit, as is shown in Fig. The course of Cj{M) in the TF 


M 


(a) 


100 


200 


300 



M 


(b) 



FIG. 2. The eigenenergy Cjk is plotted as a function of the total mass 
M for (a) y = 1, representing g > 0 and (b) p = — 1, representing 
g < 0. For the ground state (k = 0) and the radially excited states 
(k = 1, 2, 3) the course of UkiM) shows the same qualitative shape. 
In case of a fixed g > 0 and sufficiently high values of M, d)o(M) 
closely approaches the course of the TF limit which is described by 
Lo = —r/{'Ky^). In case of a fixed g < 0 a maximum value of 
the total mass is found, Maax,k (for fc = 0,1, 2, 3). The value of 
Mm!,x,k increases the higher the system is radially excited. In case 
of M < Mmax.fc two equilibrium solutions with the same number of 
nodes and for a fixed pair of parameters [g, M) can be found. 


limit is computed by using equation ( fTSl l and by converting it 
into the rescaled variables, see As predicted by ( [T8] l, the 
eigenenergy in the TF limit linearly depends on the total mass. 

For higher values of the total mass the course of Cjk{M) 
approaches a linear behavior for all fc = 0,1, 2, 3 with the 
slope of uj in the TF limit, see Fig.Flere Cjk{M) follows 
the TF limit except for a constant offset. This constant offset 
increases with the radial excitation, see Fig. 










































6 


OJfc 


10 ® 


10 ^ 


10 ”^ 



10-^ . — . — . — . ——. 

10“^ 10 10® 

M 

FIG. 3. The absolute value of the eigenenergy \uJk\ is shown as a 
function of the total mass M for g — 1,0,—1 in a log-log plot. 
Dashed colored lines correspond to g = —1, solid lines correspond 
to g = 1, and dashed black lines correspond to g = 0. Every case 
is plotted for the ground state (fc = 0) and the first three radially 
excited states (k = 1, 2, 3). For a small value of the total mass all 
graphs follow the relation oc and the non-interacting limit 
{g — 0) turns out to be a sufficient description of the system. For 
a rising value of M all graphs with a fixed g > 0 slowly reach a 
linear dependence on M with the same slope. Here follows 

the course of the TF limit except for a constant offset. 


For a fixed value g < 0 the sets of solutions of the ground 
state and the first three radially excited states possess a max¬ 
imum value of the total mass, For chosen values of 

the total mass M with M > k no equilibrium solution 
with the respective excitation can be found, see Figs.j^an d0 
Furthermore, we find Mmax.fe < .Wmax,(/c-i-i)- Thus, for values 
of M > Mmax.fc of a specific set of solutions (which is char¬ 
acterized by a fixed value of g and a fixed node number k) it 
is only possible to find equilibrium solutions (with the same 
value for g) with an increased number of nodes. For every 
total mass M smaller than Mmnx,k two solutions with k nodes 
for the same value of g can be found. Thus, for the set of 
solutions with k nodes and a fixed g we obtain two branches 
of solutions for M < Mmax,k- The branch of solutions with 
higher values of the eigenenergy uik is considered to be stable, 
the other branch contains solutions which are treated unstable. 
This issue is explained in more detail in the following section. 

For sufficiently small values of M the eigenenergy ujk of 
the unstable branch turns to be proportional to see 

Fig.0. 


FIG. 4. The eigenenergy a)*, is plotted as a function of the self¬ 
interaction factor g for M = 1. For the ground state (k = 0) and the 
first three radially excited states (fc = 1, 2, 3) Cok{g) shows the same 
qualitative shape. It exists a minimal value of the self-interaction fac¬ 
tor gmm,k < 0 for each set of solutions with k nodes. The value of 
5min,fc decreases for larger k. For a fixed gmin.fe < p < 0 it is possi¬ 
ble to compute two solutions with the same number of nodes and the 
same total mass M. 


B. The self-inter action parameter - eigenenergy relation 

In Fig. 1^ the eigenenergy ujk is plotted as a function of the 
self-interaction parameter g for a fixed value of the total mass 
M. The function {g) shows the same qualitative shape for 
the ground state (fc = 0) as well as for the first three radially 
excited states {k = 1,2,3). There exists a minimum value 
5min,fc < 0 of the self-interaction parameter for sets of solu¬ 
tions with k nodes and a fixed value of M. The value gmin,k is 
smaller the higher the system is radially excited. For a fixed 
5min,fc < 5 < 0 it is possible to compute two solutions with 
the same total mass M and the same number of nodes k', if 
k' > k. 


C. Thomas-Fermi limit of high total masses 

For a positive value of g and a sufficiently large total mass 
M the behavior of uJk{M) is well approximated by its TF 
limit. For analyzing this regime of sets of solutions in case 
of ground state solutions and their first three radial excitations 
we compute the particle densities /fe(r)^ and the appropriate 
gravitational potentials for M — 10^. Furthermore, we cal¬ 
culate /(r)^ of the corresponding TF-limit by using equation 
and by converting it into the rescaled variables ( |T^ (see 
Fig.|5}. 

While the particle density profile of the ground state is well 
described with the particle density profile of the TF limit, the 
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FIG. 5. (a) Particle density profile fk{f)^ and (b) the appropriate gravitational potentials (j)k{r) of the ground state (k = 0) and the first 
three radially excited states (fe = 1,2,3) for g = 1 and for a high total mass (M = 10^). f(r)^ in the TF limit satisfies f{r)^ = 
-kM/{ gf) sin(r/v^) and is plotted in black. The profile nearly coincides with /o(fi)^. for A: = 1, 2, 3 seem to approach /(f)^ of the 

TF limit as well, disturbed though, by gaps within their courses, which occur due to number of nodes of the wave function solutions. 


profiles of the radially excited states seem to approach that of 
the TF limit as well except for small gaps. These gaps occur 
owing to the nodes of the wave functions. In contrast to small 
M, for large M these disturbing gaps in the particle density 
profiles appear much more abrupt, see Fig. Furthermore, 
for large M the bulges in the shape of the gravitational poten¬ 
tials nearly disappear, see Fig.[^ 

Comparable results for the behavior of excited states of a 
EEC, but by assuming the EEC to be located in a harmonic 
trap, are computed in ll35l for a large self-interaction. Note 
that for the used scaling ( [T^ solutions of the GPN system at 
large total masses can be calculated from solutions for large 
self-interaction factors. This is possible, since one can calcu¬ 
late solutions for different parameters ga,Ma and gb,Mb from 
each other as long as the equation gaM^ = gtMb holds. 
Thus, solutions characterized by large values of g correspond 
to solutions characterized by large values of M. In llTSll the 
behavior of the wave functions is referred to a strong repul¬ 
sion between particles, which leads to a spatial distribution as 
even as possible. A sharper decrease of the wave function at 
the node locations allows a more even particle distribution. 

In case of the present GPN system, sharp changes of the 
wave function of radially excited states at narrow ranges 
around the nodes lead to comparably small changes at the 
other parts of the wave function. Here the kinetic term can be 
neglected and, thus, the course of the wave function again ap¬ 
proaches the TP limit. Therefore, radially excited states also 
approach the TP limit, except for a small range where gaps 
are occurring. The extension of the gaps shrinks with a rising 


total mass. 


V. STABILITY ANALYSIS 

Por the stability analysis of solutions of the dimension¬ 
less time-independent GPN system ( |20l l arguments of catas¬ 
trophe theory are applied. Likewise, this method was used 
in several publications on the stability of boson stars (e.g. 
12611221 [Ml El). A general introduction to and applications 
of catastrophe theory can be found in, e.g., If38lj40ll . Here we 
proceed similarly to the procedure employed by N. Sakai et 
al. EH. 

With the help of catastrophe theory it is possible to discuss 
the critical points of a nonlinear system. If such a system is de¬ 
scribed by a potential then its solutions are extrema of this po¬ 
tential and its critical points. In the present case we assume the 
appropriate energy functional ( fOl l is the characterizing po¬ 
tential. Solutions of the dimensionless and time-independent 
GPN system correspond to extrema of this energy functional. 

The stability of the solutions will be discussed for elements 
of the function space Dk containing the set of spherically 
symmetric functions characterized by k nodes. The stabil¬ 
ity of the ground state and the radially excited states will be 
discussed separately. 
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FIG. 6. The eigenenergy cO* is plotted as a function of fc with (a) a fixed g = — 1, (b) a fixed § — 1 and (c) a fixed M = 1. The plots indicate 
that (2)k is a strictly decreasing function of fc. 


A. Procedure to apply catastrophe theory 


As a first step we choose suitable control parameters and 
a suitable behavior variable of the system. In our case the 
behavior of the characterizing potential Z derived from the 
functional I depends on M and g. Since they can be given 
by hand, M and g act as control parameters of the system. 
The characterizing potential is described as a function of the 
behavior variable. 

To introduce a behavior variable x for our system we use 
a one parameter family of perturbed functions fk,x{f) near 
the equilibrium solution fk{f). For these functions the en¬ 
ergy functional I ( [T3] l can then be regarded as a function 
Z{x) := I[fk,x] of X. This allows us to introduce the critical 
points of the characterizing potential I{x). They are defined 
as points Xc, at which the derivation dl{x)/dx vanishes. In 
the case that the functional / is calculated with fk,x{f), then 
the equation dl[fk,x]/dx = {SI/Sfk,x)dfk,x/dx holds. Then 
we have 


SI[fk,x] 

Sfk,x 


dT{x) 

dx 


(24) 


From this it follows that 51 jSfk^x = 0 is satisfied if x = Xc- 
Then fk,xc is exactly the equilibrium solution fk- Therefore 
critical points Xc of I{x) represent solutions of the GPN sys¬ 
tem. 

It is suitable to choose /c = /(O) = x as the behavior vari¬ 
able, since it describes the system uniquely on Dk for varying 
values of the control parameters g and M. This is meant in 
the sense that one finds for every value of fc exactly one so¬ 
lution for the wave function on the function space Dk of the 
system. The corresponding family of perturbed functions can 
be specified by the relation 


M = 



In the present case it is possible to use the eigenenergy w as 
behavior variable instead of The GPN system has already 
been studied regarding Cj in the previous paragraphs. It also 


turns out that dikifc) is a strictly decreasing function of fc for 
all k (see Fig. |^. Thus, the behavior of the system does not 
change by changing its behavior variable from fc to w in any 
case. 

Now the stability analysis procedure will be introduced. 
For each state one has to; 


1. Compute the equilibrium space Aik = {dJk, M, g}, 
that contains all critical points of I[fk,uj], fk,Lc G Dk- 
Each point in Adj. represents a configuration of the self- 
gravitating EEC. 


2. Determine the degenerate points or turning points of 
the potential, which satisfy = 0, = 0. 

r j dCilk ’ OLOk 

These points are centers of catastrophes and stabil¬ 
ity changes occur here. The set of all degener¬ 
ate points in Aik is called catastrophe map = 


— =0> ii:=0}. 


dC)k 


3. Compute the energy functional I for equilibrium solu¬ 
tions fk in the neighborhood around a degenerate point 
p G E;; in order to be able to assign stability properties 
to each point m AAk- 


4. Draw the bifurcation set fk in the control space C = 
{M,g} and identify the regions of stability and insta¬ 
bility, respectively, fk is a projection of E^ onto the 
control space. 


B. Domain of stability of solutions 

Fig. 0 shows the equilibrium space AA^. This figure is a 
combination of Figs. andfor k = 3 where a qualitatively 
equal shape of the plotted functions for all k was found. Due 
to this, the shape of the respective equilibrium spaces Aik 
do not qualitatively differ, too, and we exemplarily show the 
equilibrium space Ais only. Each point contained in Ads rep¬ 
resents a solution with node number k = 3. The catastrophe 
map Efe is plotted as well. It divides the equilibrium space 
surface in an upper part and a lower part. 
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FIG. 7. The equilibrium space AI 3 contains all critical points (M , g, wa) belonging to solutions of the third radially excited state of the GPN 
system. In order to enable a better overview, tba is highlighted as a function of g for a fixed total mass M (green) and as a function of M for 
a fixed self-interaction factor g < 0 (red), g — 0 (orange) and 5 > 0 (cyan). The catastrophe set E 3 and the bifurcation set ^3 are plotted 
in black. ^3 is a projection of S3 onto the control plane C. Together with the (g = 0)-axis it divides C in three regions: S, SU and N. For 
a fixed pair of control parameters (M, g) € C one finds one stable solution, one stable and one unstable solution or no equilibrium solution, 
respectively, of the third radially excited state. 


In Fig.|^the values of energy functional I of the equilibrium 
solutions is plotted as a function of M around a point in for 
a fixed g. For each state one finds a cusp with two branches. In 
this case all points p with p G can be understood as centers 
of the fold catastrophe Il3§ti40ll . The behavior of I [fk^cj] as a 
function of uj close to p G E/c can be described by 

fb{x) = + bx, 

where x corresponds to the behavior variable uj and b corre¬ 
sponds to the control parameter M. For b = 0, fb{x) de¬ 
scribes the behavior of / at point p. For b < 0, fb{x) has 
a maximum and a minimum, which can be identified as an 
unstable and a stable equilibrium solution respectively. For 
b > 0, fb{x) has no real solutions. Regarding that, all points in 
the upper branches in Fig. [^correspond to unstable solutions 
of the GPN system and all points in the lower branches corre¬ 
spond to stable solutions of the GPN system. Since points in 
the lower branches belong to points in the upper part of Alfc 
in Fig. I^and vice versa, all points in the upper part of Alfc 
represent stable solutions of the GPN system and all points in 
the lower part represent unstable solutions respectively. 

In Fig.j^the bifurcation set ^3 is plotted in the control space 
C = {M,g}. ^3 and the M-axis {g = 0) divide C in three 


regions. For a fixed pair of control parameters (M, g) one 
finds in region S one stable solution, in region SU one stable 
and one unstable solution and in region N no solution for the 
third radially excited state of the GPN system. 


The bifurcation sets of the ground state (fc = 0) and the 
first three radially excited states (fc = 1,2,3) are plotted in 
C in Fig. [^ ^k can be calculated analytically for all states by 
using the scaling factor A, as long as for each state the value 
of one point in ^k is known, ^k is described by the equation: 
M^ax,k{9) = y/pfc/gMtnax.fe, where {M^ax,k,gk) is a known 
point in ^k- In table [n| the value of gk for Mmix,k = 1 is 
given for each state. As shown in Fig. [^ there are points in C 
for which only solutions of higher radially excited states can 
be found (e.g. in the region between ^3 of the third radially 
excited state and ^ 2 . 1,0 of the second or first radially excited 
state or ground state). The radially excited state with the low¬ 
est number of nodes has the smallest value of the energy func¬ 
tional, which also leads the smallest total energy F^tot- 
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FIG. 8. The values of the energy functional I of the equilibrium 
solutions is shown as a function of M for g = —1. Here the lower 
branches for k — 0,1,2,3 are considered to correspond to stable 
solutions of the system; the upper branches correspond to unstable 
solutions of the system. For all states the course of I shows the same 
qualitative shape. 


^^max ,k 



FIG. 9. Bifurcation set of the ground state (fc = 0) and the first 
three radially excited states (fc = 1, 2, 3) in the control space C. For 
all k shows the same qualitative shape. 


VI. CONCLUSION 

We discussed the ground state and the first three radially 
excited states of the GPN system with respect to the influence 
of two external parameters, the total mass M and the self¬ 
interaction factor g. The eigenenergy Wfe and the value of the 


state 

Qk 

ground state 

-4.1011 

1 . radially excited state 

-53.136 

2 . radially excited state 

-239.99 

3. radially excited state 

-716.18 


TABLE II. Minimum value of gk in case of M = 1 for the ground 
state and the first three radially excited states 


energy functional / are considered as functions of M or p and 
they show the same qualitative behavior for the first three radi¬ 
ally excited states and the ground state. For large total masses 
not only the ground state but also the radially excited states 
approach the TF limit. In the latter case the shape of the en¬ 
velope of the squared radially excited wave function solutions 
approaches that of the TF limit. 

By using arguments of the catastrophe theory the stability 
of the ground state and the radially excited states is discussed. 
For a fixed pair of external parameters one can find either one 
stable solution (g > 0), one stable and one unstable solu¬ 
tion {g < 0,M < Mmax.fc), or no equilibrium solution at all 
(g < 0,M > Mmax,k) for the ground state and the first three 
radially excited states. However, considering the results in 
it is highly likely that the radially excited states, which 
are considered stable in terms of the catastrophe theory, are 
intrinsically unstable for g > 0 and g < 0, M < Mmax.o- This 
was shown for a few different values of g by F.S. Guzman et 
al. by studying the time evolution of some radially excited 
equilibrium solutions, while allowing the flow of particles out 
of the numerical domain 021 . 

For g < Q and M > Mmaxfi there exists no ground state 
solution and only solutions of the radially excited states can 
be found. This leaves the fc-th radially excited state as the 
equilibrium solution with the lowest total energy if 

.^max,fc ^ ^ .^max,fc —1 • (26) 

In O 2 I there is nothing said about the intrinsic stability of the 
excited states in this case. One might expect stability for the 
lowest-energy equilibrium solution, since the total mass M 
is conserved. However, as discussed for a general relativistic 
treatment of self-gravitating BECs 1^ . this might not be the 
case. Since the lowest energy equilibrium state possesses not 
the lowest energy of all possible configurations a numerical 
discussion of these solutions in ESi revealed them to be dy¬ 
namically unstable. A similar discussion regarding dynamical 
stability for the GPN system remains to be done. 
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